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The real-time dynamics of local occupation numbers in a Hubbard model on a 6 x 6 square lat- 
tice is studied by means of the non-equilibrium generalization of the cluster-perturbation theory. 
The cluster approach is adapted to studies of two-dimensional lattice systems by using concepts 
of multiple-scattering theory and a component decomposition of the non-equilibrium Green's func- 
tion on the Keldysh-Matsubara contour. We consider "classical" initial states formed as tensor 
products of states on 2 x 2-plaquettes and trace the effects of the inter-plaquette hopping in the 
final-state dynamics. Two different initially excited states are considered on an individual plaquette, 
a fully polarized staggered spin state (Neel) and a fully polarized charge-density wave (CDW). The 
final-state dynamics is constrained by a dynamical symmetry, i.e. the time-evolution operator and 
certain observables are invariant under an anti-unitary transformation composed of time reversal, 
an asymmetric particle-hole, and a staggered sign transformation. We find an interesting interrela- 
tion between this dynamical symmetry and the separation of energy and time scales: In case of a 
global excitation with all plaquettes excited, the initial Neel and the initial CDW states are linked 
by the transformation. This prevents an efficient relaxation of the CDW state on the short time 
scale governing the dynamics of charge degrees of freedom. Contrary, the CDW state is found to 
relax much faster than the Neel state in case of a local excitation on a single plaquette where the 
symmetry relation between the two states is broken by the coupling to the environment. 

PACS numbers: 67.85.Lm, 71.10.Fd, 71.27.+a 



I. INTRODUCTION 

A unitary or anti-unitary transformation U of the 
observables and of the states of a quantum system 
leaves measur abl e quantities, such as expectation val- 
ues, invariant!^ The quantum system itself is said to 
be "symmetric" with respect to the transformation, if 
its Hamiltonian H is invariant: H' = UHU^ = H . For 
continuous transformation groups with U — exp(iA<£>) 
and parameters if, the invariance implies that the cor- 
responding generators A = A^ commute with H, i.e. 
[A, H] = 0. If the Hamiltonian is not explicitly 
time-dependent, this leads to conservation laws of the 
form ($(t)\A\^{i)) = (^(t )\A\^(t )) = const., where 

|*(i))= exp(-iff(t-i ))|*(*o)>. 

This concept might be generalized to "dynamical sym- 
metries" in the following way: Consider a transformation 
IA that leaves an observable invariant, 

A -> A 1 = UArf = A , (1) 

as well as the time-evolution operator: 

e -iH(t-t ) _^ f e -iH(t-t )y =u e -iH(t-t )y\ 

= e -iH{t-t ) _ 

This immediately implies that 

{^'{t )\e lH{t - tn) Ae- lH{t - to) \^\t )) 
= (^(t )\e lH ^- to) Ae- tH{t - ta) \^(t Q )) , 

for t > to, i.e. the time evolution is the same for different 
initial states, I*' (to)) and |*(*o))- 



Due to recent substantial progress in the field of ultra- 
cold quantum gases in optical latticespHSan almost com- 
plete real-time control of the parameters of a quantum 
system has become experimentally feasible. The study of 
ultracold fermions in optical lattices thus offers unique 
new possibilities to study non-equilibrium dynamics of 
quantum systems and the role of dynamical symmetries. 
Recently, the expansion of an initially confined quantum 
gas of fermions in the lowest band of a homogeneous op- 
tical lattice has been studied after suddenly switching off 
the confining potential, and the observed independence 
of the dynamics on the sign of the (Hubbard-type) inter- 
action could be explained by a dynamical symmetryP 

Here, we study a different situation where a dynamical 
symmetry holds locally for two different initial states but 
is broken due to their coupling to the environment. We 
consider the dynamics of a system of strongly interacting 
spin- 1/2 fermions in a two-dimensional Hubbard model 
at half-filling on a square lattice with nearest-neighbor 
hopping. This has be come accessible to studies of ul- 
tracold quantum gasesP^l The local density of spin-f 
particles, n^, is our observable of interest. It has been 
demonstrated recently with bosons that the parity of the 
number of particles can be measured with single-site res- 
olution even in a strongly correlated systempMSl For the 
Hubbard model and observables n^, we will argue that 
a particular combination IA of time-reversal, an asym- 
metric particle-hole and a staggered sign transformation 
represents a dynamical symmetry. 

An interesting question is how fast an initially pre- 
pared "classical" state builds up entanglement and cor- 
relations. This can be studied by preparing the initial 
state as the ground state of the model with isolated pla- 
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quettes p, i.e. |*(i )) = ®P laquottes b(t )), while for the 
subsequent time evolution (t > to) the inter-plaquette 
hopping V is suddenly switched on at t = t$. Double- 
well systems in one dimensiorfSHU] as we \\ as plaquette 
systems in two dimensions^ have recently been studied 
experimentally. 

A non-trivial dynamics, already for an isolated pla- 
quette, is obtained by preparing the plaquette initial 
state with the help of local fields producing a classical 
Neel state with staggered spins \p(to)) = |Neel,p) or a 
charge- density wave \p(to)) — |CDW,p) and by suddenly 
switching off the field at t = to. In case of full polariza- 
tion, these states transform into each other under 11 while 
the plaquette ground state in the absence of the fields, 
\p{to)) — |GS,p), is not invariant under U. By comparing 
the time evolution of starting from the initial plaque- 
tte product state with a global Neel excitation with the 
corresponding one for a global CDW excitation, we study 
the effects of the dynamical symmetry. In case of a Neel 
or CDW excitation that is localized to a single plaquette, 
we expect a breakdown of the dynamical symmetry due 
to the inter-plaquette coupling. 

Another goal of the present paper is a methodical one: 
Studying the real-time dynamics of a strongly-correlated 
lattice-fermion model in two dimensions far from equi- 
librium is a highly non-trivial task. Standard methods 
as dynamical density-matrix renormalization groupP^O or 
quantum Monte-Carlo^ approaches cannot be used in 
this case. We therefore resort to a cluster (plaquette) 
mean-field approach that has been developed recently for 
one-dimensional and impurity-type models, namely the 
non-equilibrium variantPHBSl of the cluster-perturbation 
theory (CPT)PHHl The NE-CPT can be characterized 
as a time-dependent multiple-scattering approach. It is 
exact for the non-interacting model with Hubbard-27 = 
and in the case of disconnected clusters, i.e. V = 0. 
In the general case, the NE-CPT starts from the one- 
particle propagator of the system of isolated clusters as 
a reference point and includes the effects of scattering at 
the inter-cluster potential V perturbatively to all orders 
but neglecting certain vertex corrections. Up to now, to 
study transient dynamics, the NE-CPT has been applied 
to one-dimensional or impurity-type models only. Here, 
our goal is to extend the theory and its implementation 
to two-dimensional models which are not easily accessible 
by other techniques. 



II. NONEQUILIBRIUM 
CLUSTER-PERTURBATION THEORY 

To discuss the main idea of the non-equilibrium 
cluster-perturbation theory (NE-CPT) and the necessary 
steps for an extension to two-dimensional lattice fermion 
models, we consider the single-band Hubbard modelP*! 
Using standard notations, the Hamiltonian that governs 



the time evolution of the system for t > to is given by: 

j 

We formally allow for an explicit time dependence of H (t) 
but will restrict ourselves to H = const, for the evalua- 
tion of the theory. In the latter case, the time evolution 
operator is given by exp(— iH(t — to))- 

The whole lattice is split into disconnected finite "clus- 
ters" . The first two terms on the r.h.s. of Eq. Q repre- 
sent the intra-cluster Hamiltonian with the intra-cluster 
hopping T H and the on-site Hubbard interaction U H . 
Note that the locality of the interaction part is a neces- 
sary ingredient for the CPT. The third term represents 
the hopping interlinking the different clusters, and the 
corresponding matrix elements are denoted by V H . The 
Hilbert space associated with an isolated cluster is as- 
sumed to be sufficiently small such that all cluster phys- 
ical quantities of interest can be calculated exactly by 
numerical means. For the subsequent calculations we ad- 
dress a two-dimensional square lattice which is tiled into 
quadratic plaquettes consisting of four sites only. 

The time evolution is assumed to start at t = to from a 
pure initial state \^>(t )) which is taken to be the ground 
state of a Hamiltonian 23p2 This has basically the same 
structure as H (t) but with possibly different parameters: 

23 = B - fiN = i T fk - l*Sjk) c) a c ka 

j,k,cr 

+ u B j2 + E v fk 4 c ** ■ ( 5 ) 

3 3,k,<? 

The chemical potential n can be adjusted to get the de- 
sired total particle number in the initial ground state. 
As the ground state of a correlated many-body system, 
the initial state must be treated by approximate means. 
Hence, to prepare for the CPT, we consider the same 
tiling of the lattice and the same corresponding split- 
ting of the Hamiltonian. Note that a generalization to fi- 
nite temperatures and thus a mixed, e.g. grand-canonical, 
equilibrium initial state could easily be described but will 
not be considered here. 

The time evolution of our basic quantity of inter- 
est, namely the expectation value (fij^(t)) of the lo- 
cal spin-f occupation number at site j, but actu- 
ally the time evolution of all one-particle observables 
A(t) = J2j,k,a a 3ka(t)cl- a c k(T with possibly explicitly 
time-dependent parameters, can be obtained from the 
non-equilibrium one-particle Green's function 

G%( Zl , z 2 ) ee -i (jc (V(^i)cL(^)) ) (6) 

via 

(A(t)) = -iJ2^k<r(t)Gl j (t + ,t-). (7) 
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FIG. 1: Contour C in the complex time plane, to is the time 
at which the system is prepared in the initial state. t max is 
the maximum time up to which observables are traced, ft is 
the inverse temperature. For a pure initial state |$(to)) we 
have /3 — > oo. 

Here, z\ and z 2 run over the Keld ysh-Matsubara contour 
C in the complex time plan c 30 l 31 l which is shown in Fig. 
[T] It consists of three branches: the upper and the lower 
Keldysh branch along the real axis (K, + and /C~, respec- 
tively) as well as the Matsubara branch M. parallel to the 
imaginary one. Here, t stands for physical (real) times 
while to — vr denotes a time point on and t £ K. . 
Furthermore, 7c is the time-ordering operator on C. Op- 
erators with a hat are given in their Heisenberg picture. 
Expectation values as well as Heisenberg time evolution 
refer to "fully interacting" Hamiltonians B and H(t), re- 
spectively. 

As opposed to a direct calculation, the main advan- 
tage to compute one-particle quantities (Eq. §f§) from 
the double-time propagator G (Eq. ^) on the Keldysh- 
Matsubara contour is that this is accessible to standard 
perturbative techniques, see Ref. [301 Here, we expand 
the full G in powers of the inter-plaquette hopping V. 
The (V" = 0)-propagator G' must be computed for the 
fully interacting model but for isolated plaquettes with 
small Hilbert space only: G"J fc (zi, z 2 ) = for lattice sites 
j and k belonging to different plaquettes. An according 
numerically exact Krylov-space technique for an efficient 
computation of G' has been introduced in Ref. 

Starting from G' , the propagator of the whole system 
can be obtained by summing inter-plaquette-scattering 
diagrams up to infinite order in V: 

G (CPT) = Qf + Qf .y.Qf + ... (8) 

For U B = U H (t) = this procedure represents an exact 
time-dependent multiple-scattering approach as all scat- 
tering paths are summed over. Furthermore, it is trivially 
exact in the case V = 0. For finite interactions and finite 
inter-cluster hopping, however, the NE-CPT propagator 
Eq. (JsJ is approximate as certain vertex corrections are 
neglected (see discussion in Ref. 21). Alternatively, the 
NE-CPT may be seen as an exact approach on the length 
scale of an individual cluster but as mean-ficld-like be- 
yond. 

A main disadvantage of the approach consists in the 
lack of any self-consistency. This is opposed to more 



elaborate cluster mean-field approaches^ as the cellu- 
lar dynamical me an-fie ld theory 3 ^ or the dynamical clus- 
ter approximation^ 3 ^ which, on the other hand, are not 
easily applicable to the non-equilibrium case. Lack of 
self-consistency or lack of variational optimization also 
implies that the approach cannot be expected to re- 
spect macroscopic conservation laws, i.e. conservation 
of the total particle number and the total spin, that 
result from the_U(l) and the SU(2) symmetry of the 
HamiltonianP^3 We therefore restrict the application 
of the method to the case of the Hubbard model at half- 
filling and vanishing z-component of the total spin, i.e. 

1 SjG P (^>(i)) = 0.5 for each plaquette p with L p = 4 
sites. In this case particle-number and spin conservation 
is enforced by manifest particle-hole and spin-inversion 
symmetry. 

A re-summation of the terms on the r.h.s. of Eq. ^ 
yields the NE-CPT equation 

G (CPT) =G , + G , -y-G (CPT) (g) 

which represents an obvious generalization of the static 
CPT25H2H to the time-dependent case. Here, all quan- 
tities are to be understood as having one spin in- 
dex, two spatial indices as well as two time arguments. 
Furthermore, C = A • B is short for CJ k (zx,Z2) = 
J2i J c dz 3 A J ; (zi, z%)Bf k (z3, z 2 )- For the inter-cluster hop- 
ping we have Vf k (zi,z 2 ) = Vjk(z{)5c(zi, z 2 ) where 5c 
is the contour delta function, and Vjk(z) = V^(t) for 
ze^ while Vjk(z) = Vf k for z e M. Note that there 
are two independent CPT equations for each of the two 
spin species. 



III. ITERATIVE COUPLING OF PLAQUETTES 

In previous studies of transient dynamics j 21 * 23 ! the NE- 
CPT has been applied to impurity-type models or small 
one-dimensional systems only. Furthermore, studies have 
so far been restricted to problems where the CPT equa- 
tion (Eq. (|9|) had to be solved only once for each parame- 
ter set. For the intended application to a two-dimensional 
lattice, however, several plaquettes must be coupled, see 
Fig. [2] for an example of an 8 x 8 square lattice with open 
boundaries. In principle this could be achieved via Eq. 
([9]) in a single step by treating the inter-plaquette hop- 
ping between all plaquettes as the perturbation V simul- 
taneously. Alternatively, however, the Green's function 
of the entire lattice may be obtained in several basic steps 
where in each step only two (possibly different) clusters 
(consisting of one or more plaquettes) are coupled at a 
time. This second procedure just corresponds to an iter- 
ation of the elementary CPT concept and also works for 
the non-equilibrium case. It is easily seen to be equiva- 
lent with the standard (NE-)CPT. However, it is clearly 
more flexible conceptually as one can make use of known 
concepts of standard multiple-scattering theory: 
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FIG. 2: (Color online) Decomposition of the Hubbard model 
on the two-dimensional square lattice into disconnected pla- 
quettes. T denotes the intra- and V the inter-plaquette hop- 
ping. U is the on-site Hubbard interaction. See text for dis- 
cussion. 



Consider the case displayed in Fig. [2] as an example and 
assume that U H (t) = U B = const, and V H \t) = V B = 
const, while Tjl(t) = T B k = const, for all sites j, k except 
for the sites belonging to the plaquette p = 1 in the up- 
per left corner. For this example, it is clearly beneficial to 
couple the cluster "1" consisting of plaquette p = 1 only, 
which is described by an intra-cluster Green's function 
Q[, with cluster "2" consisting of plaquettes p — 2, 
p = 16, which is described by the intra-cluster Green's 
function Q' 2 . Namely, Q[ \s easily obtained as the corre- 
sponding Hilbert space is small while the computation of 
Q' 2 is strongly simplified as this is an equilibrium problem. 
Q' 2 can thus be obtained by equilibrium CPT, coupling 
plaquettes p = 2, p = 16 step by step. 

As a second example let us again consider the sys- 
tem displayed in Fig. [2] but now we assume that the 
(time-dependent) inter-cluster hopping between all pla- 
quettes is the same and that all plaquettes have identical 
(time-dependent) intra-cluster parameters. In this case, 
a "cluster-doubling method" is helpful: In a first step 
plaquettes p = 1 and p = 2 are coupled by means of the 
NE-CPT. The Green's function of the resulting cluster 
"1+2" is the same as the one of cluster "3+4". Hence 
the latter is trivially obtained as a copy. In a second step, 
the clusters "1+2" and "3+4" are coupled to "1+2+3+4" 
which is again equivalent with "5+6+7+8" . This method 
can be iterated and is highly efficient. 

We may therefore concentrate on the basic step of cou- 
pling two (possibly different) clusters: Arranging the site 
indices in two blocks referring to the two clusters, we 
have 



itself while the CPT Green's function reads 

g (cpt) = ( Gn Gi2\ (n) 

\y2i V22) 

The CPT equation (Eq. @) acquires a 2 x 2 block struc- 
ture and can be rewritten in the following way: 

Gn = G[ + G[»V»G' 2 »V^ mGu, 

G22 = G 2 + G 2 »V* •G[»V»G 2 2, 
G12 = Gi • V • £22 , 

g 21 = g 2 .v^.Gn. (12) 

After solving the former two equations for the diagonal 
elements of G (CPT) , the non-diagonal ones can be eas- 
ily obtained by evaluating the latter two. For a one- 
dimensional tight-binding system the matrix V exhibits 
exactly one non-zero element, and the sums over spatial 
indices, implicit in the • operation, reduce to just a single 
term. In higher dimensions typically more than a single 
inter-cluster hopping term connect the same two clusters. 
Still the CPT equations represent a strongly sparse linear 
system of equations. Exploiting this sparsity is inevitable 
for an efficient implementation of the theory. 

IV. MODEL AND NUMERICAL APPROACH 

As an example for the application of the non- 
equilibrium CPT to a two-dimensional system and for 
a discussion of dynamical symmetries, we consider the 
two-dimensional Hubbard model on a square lattice at 
half-filling with hopping TR = TK = T between near- 
est neighbors j, k only and with U B = U H = U. The 
energy and the time scales are set by T = — 1. The 
time-evolution operator is invariant (see Eq. ^) under 
an anti-unitary transformation 

U = U X U 2 U Z . (13) 

Here, U\ is the particle-hole transformation for the spin-J, 
particles, 

Uic jX u[ = , Uic^Ul = c jX , (14) 

IA2 is a sign transformation (or (jr, 7r)-momentum boosP) 
applied to the spin-f particles, 

U 2 CjfU 2 = ZjCtf , U 2 c\^U 2 = ZjC^ , (15) 

where the sign factor Zj = ±1 for the two different sub- 
lattices of the bipartite square lattice. One easily verifies 
that UiU 2 HU\u\ =-H + const, (if Tjj - 0, Tf k = T*, 
V*l = V]^ ) where the additive constant results in an ir- 
relevant global phase in the time-evolution operator. Fi- 
nally, U3 is the anti-unitary time reversal. The Hubbard 
Hamiltonian itself is time-reversal invariant but 



for the V = propagator and for the perturbation V 



U3e -iH(t-t ) u i = e iH(t-t ) 



(16) 
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FIG. 3: (Color online) Pictorial representation of the Hamil- 
tonians B and H. The system is prepared as the ground 
state of B in the presence of a staggered field h coupling to 
the z-component of the local spin or to the charge density. 
The evolution of the respective final states is given by the 
Hamiltonian H with inter-plaquette hopping V = T (nearest- 
neighbor hopping T = —1). Here, the decomposition into 
plaquettes is artificial and defines the NE-CPT scheme. Lat- 
tice sites as well as plaquettes are numbered. Interactions and 
hoppings are indicated exemplarily only. 



Together with the sign change of H under U1U2, this 
proves the invariance of the time-evolution operator un- 
der the combined transformation U. 

We are interested in the site-dependent observable 



cLcj-f-. Obviously, 



(17) 



and thus the dynamical symmetry, Eq. holds. 

We compare the time evolutions for two different ini- 
tial states which are defined as the ground states of the 
Hubbard model B (see Eq. ([5])) where (i) the hopping be- 
tween the different plaquettes is switched off (V s = 0) 
and where (ii) an external field is switched on. The "spin- 
excitation state" |Neel,p) is the plaquette ground state 
in the presence of a staggered magnetic field: 



B B — h Zj (n^ 



(18) 



where Zj = (— 1) J . The "charge-excitation state" 
|CDW,p) is obtained with the help of a staggered po- 
tential: 



B^B — h } j Zj (n 

3 



(19) 



The final-state dynamics for t > to is governed by the 
Hubbard Hamiltonian H with the hopping between the 
different plaquettes switched on, i.e. — V = T 

between nearest neighbors j and k, but with the field 
switched off, i.e. h = 0. 

The actual NE-CPT calculations are performed for a 
finite system with 6x6 sites and open boundary condi- 
tions. For this system size the numerical computations 
can be performed conveniently on a standard desktop 
computer. While the initial states are taken as simple 



product states (Fig. [3j left) and are easily computed by 
means of exact diagonalization, an approximation must 
be used to determine the final-state dynamics: The inter- 
plaquette hopping V — T is treated perturbatively to all 
orders within the CPT, see Fig. [3j right. One plaquette 
at a time, in the order indicated by the numbers, is cou- 
pled in a sequence of 8 NE-CPT steps. We have checked 
that the results do not change within numerical accuracy 
when using a coupling scheme with a different ordering. 

As V B = for the initial state, there are no potential- 
scattering vertices on the Matsubara branch of the con- 
tour. Consequently, only the two Keldysh branches must 
be taken into account in the CPT equations, Eqs. (12 1. 



For the numerical solution of the corresponding integral 
equation, we introduce a time discretization with a finite 
time step At. Eqs. ( 12 ) are considered as linear systems 



of equations for the different components of the CPT 
propagators. This component decomposition is advanta- 
geous as it allows an efficient use of standard quadrature 
rules. Details are discussed in the appendix [] Apply- 
ing the trapezoidal rule, a time step At = 0.05/|T| has 
turned out as sufficient for converged results. 



V. NUMERICAL RESULTS AND DISCUSSION 

We first consider a global excitation, i.e. h is finite 
on all sites for the initial state, and the j-sums in Eqs. 
(18) and (19) extend over the entire lattice. In case of 
a small field strength h, we find clearly different time 
evolutions of the local spin-f occupation number (njf(t)) 
for the two different initial states |Neel) = <g) p |Neel,p) 




FIG. 4: (Color online) Time evolution of the local spin-f 
occupation numbers at the sites j = 1, j — 4 for the 
system displayed in Fig. [3] (with the staggered field applied 
to every plaquette). Inter-cluster hopping V = T, Hubbard 
interaction U = 8\T\. Energy and time scales are set by 
T = — 1. Lines: initial "spin-excitation state" |Neel) obtained 
with h = 100 in Eq. (l8i. 
state" |CDW), h = 100, 
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see Eq. (l9l). 
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FIG. 5: (Color online) Time evolution of the local spin-f occupation numbers at the sites j = 1, j = 4 as in Fig. [4] 
(U — 8\T\, T = —1) but with the initial field (strength h = 100) applied to the sites in plaquette p — 1 only. The initial states 
of plaquettes p = 2, p — 9 are the respective ground states |GS,p). Top (a,b): Intra-plaquette dynamics only (V = 0). 
Bottom (c,d): Plaquettes coupled via NE-CPT (V = T). Left (a,c): initial "spin-excitation state" on p = 1, i.e. |Neel, 1). Right 
(b,d): initial "charge-excitation state" on p = 1, i.e. |CDW, 1). 



and |CDW) = ® p |CDW,p). In fact, the dynamical 
symmetry is expected to hold for h — > 00 only. In 
this limit, a fully polarized Neel state | \., t>4jt) an d 
a CDW state |0,f4, 0, tl) are created on each plaque- 
tte, respectively, which are mapped onto each other: 

U,t4,t)=w|o,t+,o,tO. 

The results for a large field ft, = 100 are shown in Fig. 
[i| At t = t a = we find (n jt (t)) « for j = 1 and j = 3 
and (njf(t)) ~ 1 for j = 2 and j = 4 (see Fig. [3] for the 
labeling of the sites) for both, the spin- and the charge- 
excitation state. Note that the occupation numbers of 
the spin- J, particles are fixed by particle-hole and spin- 
inversion symmetry, i.e. (nj±(t)) = l — (rij-f(t)} (Neel) and 
(riji(t)) = (rijf(t)) (CDW). Hence, total particle-number 
conservation and, in the same way, conservation of the 
z-component of the total spin is enforced. This sym- 
metry has also been verified numerically. Furthermore, 
spatial symmetries (see Fig. [3]) enforce equal occupations 
(fij-f(t)) for the sites j = 2 and j = 4 while sites j = 1 
and j ' = 3 are inequivalent. As can be seen from Fig. [4j 
this is respected by the NE-CPT. The approximate ap- 
proach is also seen to respect the dynamical symmetry: 
The results obtained for the two different initial states 
I Neel) and |CDW), displayed as lines and dots in Fig. 
[4j respectively, are almost equal in the entire time inter- 



val studied. Remaining discrepancies are small and are 
due to residual numerical errors in the solution of the 
CPT equation as well as due to tiny deviations from full 
polarization in the initial state at h = 100. 

Fig.[5](a) and (b) displays the time evolution of (rijf(t)) 
starting from initial Neel and CDW states (h = 100) 
for the sites in the isolated plaquette p = 1: The inter- 
plaquette hopping stays switched off. Due to the dynam- 
ical symmetry, the time evolution is the same in both 
cases (a) and (b) . In Fig. [5] (c) and (d) the correspond- 
ing results are shown for an initial Neel and CDW state 
prepared on plaquette p — 1 only (h = 100) while pla- 
quettes p = 2, p = 9 are assumed to be in their 



singlet ground state |GS,p), i.e. the j-sums in Eqs. (18) 
and ( 19 ) extend over p = 1 only. Contrary to (a) and 
(b), the inter-plaquette hopping is switched on (V = T) 
for times t > 0. 

We first note that U transforms the initial states on 
plaquette p = 1 into each other, |Neel, 1) = U\CUW, 1), 
but does not leave the singlet ground state on the other 
plaquettes invariant, i.e. |GS,p) ^ U\GS,p). Conse- 
quently, the dynamical symmetry for the two different 
initial states is lost and different time evolutions are ex- 
pected for finite V. This is nicely seen in Fig. [5] by com- 
paring (c) with (d). It is striking, however, that in case 
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of the initial Neel state (c) the time evolution of (n^(t)) 
for sites j in plaquette p = 1 does not deviate much from 
the pure intra-plaquette dynamics shown in (a). On the 
other hand, starting from the initial CDW state (d), we 
find strong inter-plaquette effects dominating the dynam- 
ics almost immediately after the quench. 

For an explanation of these findings, we first concen- 
trate on the intra-plaquette dynamics. Here, it is suf- 
ficient to discuss the case of the Neel initial state only: 
As is seen in the figure, the full spin polarization quickly 
decays but at t T (U)/2 ss 9.3 for U = 8 an almost fully 
polarized state with a reversed sign of the local spin is 
recovered. Note that t T (U) is different from and actually 
much smaller than the exact recurrence time of the fi- 
nite quantum system. However, it does represent the ex- 
act recurrence time in the strong-coupling limit U — > oo 
where the low-energy sector of the plaquette Hamiltonian 
is exactly mapped onto an antiferromagnetic Heisenberg 
plaquette. 

Consider a dimer model rather than a plaquette for a 
moment: For U — > oo, the Neel initial state |Neel,p) = 
| t; I) couples to the singlet ground state and to the 
triplet excited state of the effective Hei senberg model 
H low = JS 1 S 2 , namely | t,l) = V^(\S = 0, M = 
0) + 15 = 1,M = 0)). The time evolution of (% t (*)> = 
0.5 + (Sj^ z (t)) is thus governed by a single frequency ui 
given by the excitation energy u> = E(S = 1) — E(S = 
0) = J = 4T 2 /U which translates into a recurrence time 
t T (U) = 2tt/lo = ttU/(2T 2 ). 

For the plaquette p = 1, |Neel, 1) turns out to be given 
by a linear superposition of three energy eigenstates for 
U — > oo. The resulting three excitation frequencies ^2,3, 
however, are integer multiples of the smallest one, say 
uj\. Rescaling the corresponding recurrence time t r (U) — 
2ir/uii oc U gives us t r (U = 8) « 12.6 which is close 
but not equal to the observed t T (U) ss 18.6. This is 
an effect of residual charge fluctuations which are absent 
in the pure spin model and small but non-negligible for 
[7 = 8. Charge fluctuations are also responsible for the 
small wiggles which can be seen in Fig. [5] (a) on a time 
scale 2tt/U < 1 and which are absent in the Heisenberg 
limit. 

Having understood the physical origin of the different 
structures in case of isolated plaquettes, we can interpret 
the results in Fig. [5] (c) and (d). The initial Neel state 
on p = 1 has a mean energy (Neel, l|i?|Neel, 1) ~ J and 
thus preferably couples to the low-energy eigenstates of 
the full model while for the initial CDW state we have 
(CDW, l|ff|CDW, 1) - 2U, i.e. the latter mainly couples 
to highly excited states. Note that the energy spread 
AE = yJ{H 2 ) - (H) 2 ss 2.8 is almost the same for both, 
the Neel and the CDW initial state. 

For the Neel initial state this implies that the overall 
trend of (n^(i)) is much slower as compared to the CDW 
initial state. Therefore, on the time interval displayed in 
panel (c), the dynamics still follows the intra-plaquette 
dynamics (a) which is governed by the time scale 2tt / J. 
The much less pronounced structures on the time scale 



2ir/U that are induced by residual charge fluctuations, 
on the other hand, are almost immediately affected by 
the scattering at the inter-plaquette potential V. This is 
clearly visible, e.g. by comparing the results for inequiv- 
alent sites j = 1 and j = 3 in (c). Relaxation of (jij\{t)) 
to its equilibrium value (rijf) — 0.5 is thus expected for 
times much larger than 2it / J . 

In case of the initial CDW state, the overall features of 
the intra-plaquette dynamics (b) are dominated by the 
same time scale 2n / J due to the dynamical symmetry. 
For coupled plaquettes (d) where the dynamical symme- 
try is absent, however, the dominant time scale is rather 
given by 2ir/U <§; 2ir/ J . This explains the observed much 
faster relaxation of (%f(t)) on a time scale larger than 
2-k/U but smaller than 2ir/J. 

We conclude that the largely different time evolutions 
displayed in (c) and (d) result from a clear separation of 
time scales while the dynamical symmetry enforces equal 
behavior for decoupled plaquettes as seen in (a) and (b). 
This interpretation is corroborated by calculations at a 
much weaker interaction strength U = 2 (not shown) 
where nominally J — U but where there are actually 
no well-defined local moments. Here, the relaxation of 
(rij^(t)} is indeed equally fast for both types of initial 
states. We also note that in case of a global excitation 
(see Fig. [4]), where the dynamical symmetry is intact, the 
results for (n^(i)) again closely resemble those obtained 
for the isolated plaquette, see Fig. [5] (a) and (b). 

Finally, our cluster mean-field results for the two- 
dimensional Hubbard model can qualitatively be com- 
pared with numerically exact data obtained by means 
of time-dependent density-matrix renormalization-group 
methods (DMRG) for the one-dimensional Hubbard 
modelJMMI j n j^ e f [Jgj a i oca i doublon-holon excitation 

of the ground state of the half-filled model has been con- 
sidered. The doublon is found to delocalize quickly on a 
time scale of a few 1/|T|, very similar to doublon der- 
ealization in an empty- band system,^ but, on a scale of 
several tens of 1/|T| , does not completely decay as higher- 
order many-body scattering processes are necessary in 
the strong-coupling regime to break up the repulsively 
bound fermion pair consistent with energy conservation. 
Note that within the NE-CPT there is no direct access 
to two-particle expectation values. Nevertheless, the fast 
derealization of the initial charge excitation is qualita- 
tively consistent with our results in Fig. [5] (d). 

For the homogeneously excited initial state |CDW) 
where a doublon is present at every second lattice site, 
DMRG 39 again predicts the doublons as very stable in 
the strong-coupling regime, i.e. the total double occupa- 
tion per site (D(t)) relaxes to a constant which is close to 
its maximum initial value (D(0)) = 0.5 (an exponential 
decay is only found in the presence of a finite nearest- 
neighbor Coulomb interaction). Our results (Fig. [4]) are 
again qualitatively consistent with the DMRG data: At 
t\ 9.3 we find the local occupation numbers close to 
unity (j = 1,3) or zero (j — 2,4). Here, we can easily 
estimate (D(h)) « Ej Gp <%t(*i))<%4.(*i))/£p ^ 0.5. As 



the occupation- number dynamics (rij a (t)) does not differ 
very much from the dynamics for an isolated plaquette on 
the time interval considered here, we expect the same for 
(D(t)). For the isolated plaquette, we in fact find (D(t)) 
slightly oscillating around a value smaller but close to 
0.5. 



VI. CONCLUSIONS 

The generalization of the cluster-perturbation theory 
to systems with a non-trivial real-time dynamics of an 
initial non-equilibrium state provides a conceptually ap- 
pealing approach that has been demonstrated as nu- 
merically feasible even for two-dimensional lattice mod- 
els of strongly correlated fermions. The NE-CPT rep- 
resents a non-self-consistent cluster mean-field-type ap- 
proach where the effect of the scattering at the inter- 
cluster potential on the self-energy is neglected. For 
the particle-hole and spin symmetric Hubbard model, it 
nevertheless respects macroscopic conservation laws and 
has previously been showrP^l to be able to reliably de- 
scribe relaxation dynamics on short time scales for one- 
dimensional and impurity-type models. 

The present study has addressed important numer- 
ical improvements that are necessary to study two- 
dimensional models: First, for a numerically efficient 
evaluation of the NE-CPT equations on the Keldysh- 
Matsubara contour, we have rewritten and implemented 
the NE-CPT equations for the different components 
of the non-equilibrium Green's function. This leads 
to a considerable increase of the numerical accuracy 
and the accessible times. Secondly and more impor- 
tant, however, we have demonstrated that standard con- 
cepts of multiple-scattering theory, known from ab initio 
electronic-structure theory, can successfully be applied to 
the non-equilibrium case. The iterative scheme of cou- 
pling of only two equal or different building blocks at a 
time represents an efficient and flexible concept that is 
applicable to the two-dimensional case but may also be 
useful in different and even more complicated situations, 
e.g. in cases with fewer spatial symmetries. 

As a concrete example, the dynamics of the Hubbard 
model on a 6 x 6 square array has been considered for 
initial states prepared as plaquette-product states. We 
have studied the time evolution of the spin-dependent 
local occupation numbers to understand relaxation ef- 
fects induced by inter-plaquette correlations building up 
in the course of time. Their real-time dynamics in the two 
states |Neel) = ® p |Neel,p) and |CDW) = ® p |CDW,p) 
with a global staggered spin or charge excitation is found 
to be identical. This is enforced by a dynamical sym- 
metry, i.e. an invariance of the time-evolution opera- 
tor and of the considered observables under a combined 
spin-asymmetric particle-hole, staggered sign, and time- 
reversal transformation U, which is found to be respected 
by the NE-CPT. This dynamical symmetry prevents a 
fast relaxation of the CDW initial state which was ex- 



pected to take place on a short time scale 2ir/U related 
to charge dynamics. 

If, on the other hand, the dynamical symmetry is bro- 
ken by preparing the spin or charge excitation locally on 
the plaquette p = 1 only while the rest of the system 
is given as a product of plaquette ground states |GS,p), 
a clear separation of energy and time scales is observed 
in the relaxation dynamics: The occupation numbers in 
the spin-excited plaquette (|Neel,p = 1)) show an overall 
slow coupling to the environment on the scale 2n/ J, with 
some fast but small oscillations due to residual charge 
fluctuations on top. Contrary, a charge dynamics on the 
scale 2ir/U <C 2ir/J is seen to result in a fast relaxation 
of the of occupation numbers for the charge-excited state 
|CDW,p = 1). 

One might speculate on the long-time limit that is 
not accessible to real-time Green's-function-based ap- 
proaches as the NE-CPT: For the global spin and charge 
excited states, the dynamical symmetry enforces equal 
time evolutions although the total energies per site are 
largely different (~ U/2) which therefore implies differ- 
ent thermal states characterized by different tempera- 
tures eventually. In fact, the dynamical symmetry trans- 
lates into a symmetry of the thermal states as H t— > 
—H + const, under IA implies a transformation of the 
(normalized) thermal mixed state p i— > p' with a negative 
temperature T h-> — T. For the locally excited states, on 
the other hand, the higher energy of the local charge ex- 
citation is expected to dissipate to the bulk of the system 
resulting, in the thermodynamical limit, in equal temper- 
atures and equal thermal states eventually. 
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Appendix: Component decomposition 

The • operation occurring in the CPT equation 
(Eq. ([9])) includes a convolution with respect to time vari- 
ables. The corresponding integration path is the contour 
C shown in Fig.[T] If C was considered as a whole, the er- 
ror due to the time discretization in the numerical treat- 
ment of the integrals, e.g. via a Newton-Cotes formula, 
would be dominated by the characteristic discontinuities 
of non-equilibrium Green's functions originating from the 
time ordering along C, see Eq. (|6j. Smaller errors or 
larger time discretization steps and thus a larger t max , 
however, are attainable by splitting the non-equilibrium 
Green's functions into certain components based on the 
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different branches of C. Exemplified by G and omitting 
spin and spatial indices for simplicity, one may define 



G<(t u t 2 ) = G(t+^) , 
G ,n (t 1 ,r 2 ) = G(<±,i -iT- 2 ) , 
G r "(Ti,i 2 ) = G(t -i Tl ,t±) , 



(A.la) 
(A.lb) 
(A.lc) 
(A.ld) 



G MG (n,T 2 ) = -i (c{t -in) c\t Q -iT 2 )) , (A.le) 

G ML (n,r 2 ) ee i (c*(t -iT 2 ) c{t -i n )) , (A.lf) 

G ret (h,t 2 ) ee 0(h-t 2 ) (G>(t u t 2 )-G<(h,t 2 )) , 

(A.lg) 

G adv (i 1; t 2 ) ee -Q(t 2 -tt) (G>(t u t 2 )-G < (t 1 ,t 2 )) . 

(A.lh) 

Besides the distinction between "greater" and "lesser" 
Matsubara components G MG and G ML , the component 
definitions are well-known in non-equilibrium Green's 
function literature.^ For our present purposes, G rct and 
G adv are auxiliary quantities only, all necessary informa- 
tion on G is already contained in the components G > , 
G ML . We furthermore write 

(A C1 o*J B c *)( Zl ,z 2 ) ee f U dt 5 A c Hz 1 ,t 5 )B c Ht 5l z 2 ) , 

Jt 3 

(A.2) 

and 

(A C1 B c -){ Zl ,z 2 ) = -% j ^ d^A^{z U 7t)B^{^,Z2) ■ 

Jt 3 

With this, a contour convolution C{z\,z 2 ) = 
J c d Z 3 A(zi, z 3 ) B(z 3 , z 2 ) can be split into the different 
components and reads: 



These convolution rules represent an adaption of the 
Langreth rules^l for our purposes and are similar to cor- 
responding rules described in Ref . 2U Application to the 
differential equation of motion for G yields generalized 
Kadanoff-Baym equations!^ Utilizing the component de- 
composition and also evaluating the unit step functions of 
the retarded and advanced components, we find all inte- 
grands being continuous functions of the respective time 
argument. Furthermore, causality is made explicit by the 
(renewed) integration boundaries. Although the convo- 
lution G does not behave like a proper Green's function 
in all aspects, it can play the role of A or B in a next 
convolution step, e.g. the components G c are continuous 
themselves. 

For the NE-CPT, we actually have to solve a problem 
of the form A • B = C where A and C are given. To 
this end, we rewrite the CPT equation for the different 
components and treat the resulting problem as a system 
of linear equations (1 - G' • V) • G (CPT) = G' which 
must be solved for G^ CPT ^. Because of the homogeneity 
of the Matsubara Green's function, solving Eq. (A.4e) 
with t 2 = is sufficient to get B MG and B ML . Thes e 
components in turn provide us with B n via Eq. (A4c |. 
Equipped with B adv from Eq. ( |A.4h| , Eq. ( |A.4df ~can 
be solved for B r . Both components are prerequisites 
for finally obtaining B> and B < from Eqs. (A.4a) and 
(A.4b|, respectively. The lesser component is our main 



quantity of interest as it enters Eq. The mentioned 
solution steps can be performed separately for all spa- 
tial/temporal indices k and z 2 of BJ k (zi,z 2 ) with the 
additional benefit of recurring coefficient matrices. Fur- 
thermore, use can be made of the fact that coefficient 
matrices representing A lct and A &Av are triangular. 



C> = A> tm o ax B adv + A Ict *"o x B> + A" * B r , (A.4a) 

to to 



G< = A< tm o ax B adv + A Ict * n o x B<+A"iB , (A.4b) 

to 



.1 „r 
* 




G n = A TCt *"o x B n + A" * S ML + A" * B MO , (A.4c) 



G r = A r tm o ax B adv + A MG * B r + A ML * B r , (A.4d) 

to ri 

Tl >T 2 

G MG = A MG 1 B ML + A MG * B MG + A M1 - I B MG , 



n<T 2 

C ML = A MG * B ML + A Mlj 1 B ML + A ML * B MG 

Tl T 2 



QtBt = A mt Y B ret , 
to 



G adv = A adv o ax B 
t 



adv 



(A.4e) 

(A.4f) 
(A.4g) 

(A.4h) 



Here, we have also omitted time variables except for those 
appearing as integration boundaries. r a denotes the a th 
argument of the component C' c . 
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